Convection in multiphase flows using Lattice Boltzmann methods 
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We present high resolution numerical simulations of convection in multiphase flows (boiling) using 
a novel algorithm based on a Lattice Boltzmann method. We first validate the thermodynamical and 
kinematical properties of the algorithm. Then, we perform a series of 3d numerical simulations at 
changing the mean properties in the phase diagram and compare convection with and without phase 
coexistence at Ra ^ 10^. We show that in presence of nucleating bubbles non-Oberbeck Boussinesq 
effects develops, mean temperature profile becomes asymmetric, heat-transfer and heat- transfer 
fluctuations are enhanced. We also show that small-scale properties of velocity and temperature 
fields are strongly affected by the presence of buoyant bubble leading to high non-Gaussian profiles 
in the bulk. 

PACS numbers: 47.20.Bp,47.55.-t 



Thermal convection, the state of a fluid heated from 
below and cooled from above, is an ubiquitous phenom- 
ena in nature, present in many industrial and geophysi- 
cal applications both at micro and macro-scales [1]. It 
is also challenging from the theoretical point of view, 
due to its extremely reach and different regimes rang- 
ing from intricate pattern formations at small tempera- 
ture difference between bottom and top plates (i.e. mod- 
erate Rayleigh number) to extremely turbulent behav- 
ior where heat transfer and its adimensional definition 
(i.e. Nusselt number) is dominated by bulk or bound- 
ary layer physics (or by both, see e.g. recent reviews 
[2]). Most of the time thermal convection is studied in 
its simplest version, the so-called Oberbeck-Boussinesq 
(OB) approximation, where a single phase -unstratified- 
fluid is present with constant material properties. Com- 
pressibility is also neglected except for buoyancy forces. 
Needless to say, in many situation some, or all, of the 
above assumption breaks down and one enters in the 
realm of Non-Oberbeck-Boussinesq (NOB) convection. 
Deviations from OB can arises in many different way. 
Two notable cases are (i) the presence of stratification, 
as in many geophysical applications (ii) and under boil- 
ing conditions, i.e. when the parameters excursion inside 
the convective cell allows for phase coexistence [4] [5] . 

In this paper we address the latter case. We study ther- 
mal convection in a 3d cell in a high turbulent regime 
where large bubbles (larger than the turbulent viscous 
scale) can nucleate in the layer close to the bottom wall 
with a non-negligible heat-exchange between liquid and 
vapor. To do that, we present, validate and apply a 
novel numerical scheme based on a diffuse interface Lat- 
tice Boltzmann method (LBM) O [?]• We are not re- 
stricted to treat bubbles as point-like [8] and we fully 
resolve the thermo hydrodynamical properties of the gas 
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FIG. 1: Check of Clausius-Clapeyron relation dP/dT — 
As/Av in our numerical scheme at varying T/Tc. P is the 
equilibrium pressure at coexistence temperature T, v = 1/p is 
the specific volume, Tc is the critical temperature and s(T, p) 
is the specific entropy. Bottom inset: latent heat, A = TAs vs 
T. Top inset: Bubbles are in blue. Regions with high temper- 
ature are in red. The system has no-slip velocity at bottom 
and top walls and it is periodic on the horizontal directions. 



and liquid phases. Beside the methodological aspects, 
we also address physical questions connected to the en- 
hancement/depletion of heat flux in presence of bubbles, 
statistics of mean global properties as well as small-scales 
effects for both velocity and temperature fluctuations. 
We present two series of high-resolution numerical sim- 
ulations up to 512^ collocation points at Ra ~ 10^ with 
and without phase coexistence, such as to be able to di- 
rectly compare on the same geometrical set up the effect 
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FIG. 2: Top panel. Phase space T — p equilibrium curves 
(solid lines) superposed with a scatter plot of temperature 
and density values measured under boiling condition (both 
made dimensionless using the critical values, Tc and pc). No- 
tice the presence of bubbles with different temperature inside 
the volume. Different symbols corresponds to measurements 
taken in top, bottom or bulk region. The horizontal dashed 
line correspond to top, Tu^ bottom, T^, and m ean, T m, tem- 
peratures. Inset: mean temperature profile, T[z)/Tc vs z 
(in lattice units) for boiling and non-boiling conditions. The 
straight line corresponds to the adiabatic slope. 
Bottom panel. Bulk contribution to the heat flux normal- 
ized to its diffusive value, pUuz/[i^{Td — Tu)/L\^ (Nusselt) for 
boiling and non-boiling case at comparable Rayleigh number. 
Inset: pdf of pUuz normalized to have mean area and mean 
variance for both boiling and non-boiling cases. 



density contrast, interface thickness, viscosity and ther- 
mal diffusivity to realistic situations. This is why the 
interplay between numerical investigations and experi- 
ments is crucial in this field. 

The equations of motion describing a non-ideal fluid in 
presence of thermal fluctuations are: 

dtpui^dj{puiUj) = -diP^dj{ii{diUj^djUi))^gpz (1) 

where ja = pu is the molecular viscosity, g is the gravity, p 
is the local fluid density and P(p, T) = Po{p, T)-\-Pni{p) 
is the non-ideal pressure. Pressure is fixed by the equa- 
tion of state and it is made of two terms, the ideal part 
Po(p, T) — pT and the non ideal part which in our LBM 
system reads: Pni{p) = Gexp(— 2/p) (see below). The 
equation for the internal energy, U = CyT + J dpPNi / p^ 
is given by: 



pDtU 



(2) 



where n is the thermal conductivity. The above equa- 
tion can also be rewritten as one of the two following 
equivalent expressions: 



CppDtT - pTDtP 
CypDfT ^ PodjUj = 



K.djjT 



(3) 



where Dt stands for the material derivative, Cy is the 
specific heat at constant volume, Cp and /3 = —{dTp)/p 
are the specific heat and compressibility at constant pres- 
sure, respectively. The above set of equations tends to 
the usual OB system when the fluid is considered single 
phase, incompressible and both p^ k are constant [6\ In 
table 1 we report the characteristic values for all relevant 
parameters, with and without boiling. 

Because of bubble nucleation/evaporation, a key role 
is played by the DtP term in ([3|. Take for example a 
convective cell of eight L with imposed temperature, T^^, 
at the bottom wall and Tu at the top wall. Then, the heat 
balance across an horizontal layer at distance z from the 
bottom wall is given by the expression: 



dtpU\z -^dzpUvz - i^dzT\z = 



-PdjUj\ 



(4) 



where with {•)\z we intend a spatial average at fixed z. 
In a stationary situation, we can define a z-dependent 



dimensional Nusselt number Nu{z) = pUvz 
which satisfy an integral constraint 



■ ndzTV 



of boiling on convection. With respect to experimental 
studies, numerical simulations offer the unique advan- 
tages to allow access to all quantities without affecting 
the fluid dynamics and to confine the fluid inside "ideal" 
surfaces (i.e. perfect thermal properties at the wall, per- 
fect smoothness of the boundaries etc..) On the other 
hand, a limitation consists in the difficulty to reach high 
Rayleigh numbers and to push the physical parameters as 



Nu{z) - Nu({)) 



dz' PdjUj\z' 



(5) 



With the above definition, the Nusselt number is not any- 
more constant throughout the cell, we may exchange heat 
by nucleating and evaporating bubbles or by simple com- 
pressible effects inside each phase. 

Algorithm. The numerical algorithm used is based 
on discrete kinetic models [7 . The starting point is a 
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AT Acp Xg XI ^ Rci A/3 Pr 


BOILING 
NO-BOIL. 


0.226 1.2 0.008 0.0018 0.0165 3 x 10^ 1.6 9 
0.230 0.2 — 0.0018 0.0165 2 x lO'^ 0.1 9 



TABLE I: AT, Acp, A/3: values of temperature, liquid heat 
capacity and liquid compressibility difference between the two 
walls (all normalized with their respective values at the center 
of the cell, Zc). Xi,g — K/{cp{zc)pi,g{zc)): thermal diffusivity 
of liquid (1) and gas (g). u: kinematic viscosity. Rayleigh Ra 
and Prandtl Pr numbers are evaluated at Zc and in the liquid 
phase: Ra = ^^(^c)^"(-^^/-^-7ad(.e)) ^ _u_ ^^^^^^ 

L = 512 is the cell height (in grid units) and 7ad = /3Tg/cp 
is the adiabatic gradient. The Jacob number quantifying the 
ration between the sensible heat and the latent heat [8] is 
Ja ^ 3. 



standard coupled mesoscopic dynamics described by \T, 

fi{x^ci,t^l)-fi{x,t)= -^{fi-f^''^){x,t) (6) 

gi{x^ci,t^l)-gi{x,t)= -^{gi-gl'^'^Xx.t) (7) 

where fi{u^t)^ giix^t) stand for the probability density 
functions to find at ix^t) a particle whose kinetic velocity 
belongs to a discrete and limited set q (with / = 1,19 
in the 1)3(519 LBM adopted here [7 ). Density, momen- 
tum and temperature are defined as coarse-grained (in 
velocity space) fields of the distribution functions 

p=Y.fi pu = Y^cifi pT = Y^gi. (8) 

I I I 

The local kinetic equilibria //^^^ {u' , p) and g\^^^ (n, F, T) 
are usually expanded in suitable polynomial basis [10] in 
such a way that a Chapman-Enskog expansion [9] leads 
to the equations for density, momentum and tempera- 
ture ([T])-(|3|: the streaming step on the left hand side of 
(|6| reproduces the inertial terms in the hydrodynamical 
equations, whereas dissipation and thermal diffusion are 
connected to the relaxation (towards equilibrium) prop- 
erties in the right hand side, with v and n related to the 
relaxation times Tj^^ [7|. Non ideal thermodynamics is 
obtained by a well controlled procedure shifting the ve- 
locity in the equilibrium distribution, = u -\- r^F/ 
with a forcing term mimicking the effect of an internal 
pseudo-potential [9l|10]. In particular, we adopt the stan- 
dard form: 

N 

F = -gY, w{\ci\^)ci^[p{x)]^[p{x + ci)] (9) 

where ^ is a parameter dictating the overall strength 
of the non-ideal interactions. The weights w{\ci\^) are 
used to enforce isotropy up to the 4th order in the ve- 
locity tensors [11]. The pseudo-potential, encom- 
passes the macroscopic effects of both long-range attrac- 
tion and short-range repulsion. Although various choices 



have been presented for the choice of ijj [p] [131 E] ? here it 
is crucial to set it to = exp( — 1/p), in such a way to 
reproduce the thermodynamic consistency on the lattice 
(see fig. [T]). Furthermore, to reproduce the correct 
ideal part of the pressure, Pq = pT, a coupling between 
fi and gi populations in ([3| is needed. To this end, we 
used a recent proposal [16 by plugging the dynamical 
temperature T extracted from the population gi in the 
equilibrium distribution of ([7|): in the limit of vanishing 
interaction, this is equivalent to impose a second order 
momentum of //^^^ equal to //^^^c^c;^ = pT5ij-\-pUiUj. 
Finally, in order to reproduce exactly the divergence 
term, P^djUj^ in (|3|, we found necessary to add a proper 
counter term to the evolution of gi populations in ([t]), 
as proposed in [15 . As a result, we ended with a LBM 
scheme able to reproduce in the hydrodynamical limit 
the NS equations ([l][3| with a non-ideal Pressure tensor 
and a consistent definition of latent heat (see fig. [T]). 

Single point quantities. In fig. Q we shown a 
scatter plot of T(x, t) vs p(x, t) for a boiling cell. As 
one can see most of the volume is at thermodynamical 
equilibrium^ superposing with the equilibrium curves 
in the T — p phase space. The presence of bubbles is 
clearly detected by the spots concentrating along the 
vapor branch and it is also interesting to notice that the 
corresponding bubble temperature is always larger than 
the mean temperature in the cell, indicating that bubbles 
are transferring temperature upwards very efficiently. 
Moreover, the temperature profile across the cell, T(z), 
becomes slightly asymmetric in presence of bubbles, 
a phenomenon also observed in other liquid-like NOB 
systems [17 . Breaking of the top-down symmetry must 
not surprise. In particular, P is not constant across the 
cell (i.e. density and temperature fluctuations are not 
strictly proportional as in OB) and cp decreases going 
from bottom to top. Both effects may have an impact 
on the averaged profiles as discussed and observed also 
in [17 . Here, the temperature mismatch between the 
values at the center and the mean temperature is 1% 
(inset of top panel in fig. |2|. Notice also, in the same 
plot, that T{z) agrees with the expected profile given by 
the adiabatic gradient, due to the presence of a small 
stratification. 

In the bulk, the heat flux ([5| is dominated by the con- 
vective term pUuz- In bottom panel of fig. ^ we com- 
pare the Nusselt number for the boiling and non-boiling 
cell at comparable Rayleigh. Two effects show up. First, 
heat flux is enhanced. Second, fluctuations around its 
mean profile are larger in presence of bubbles. We in- 
terpret this as a clear signature of the importance of the 
bubble dynamics in transporting heat between the two 
walls. This is the combined effect of temperature entrain- 
ment inside bubbles leveraged with the buoyancy effect 
that drift bubble upward much more efficiently then for 
plumes in single-phase convection. Because bubbles are 



rare in our system, this also implies an increase in heat 
flux fluctuations, as can be seen in the inset of bottom 
panel in flg. ([2| where we show the probability density 
function (pdf) of the heat flux measured only in the bulk 
cell. Clearly, the right tail are enhanced, due to bub- 
ble buoyancy. An open question is to see whether this 
enhancement of heat-flux is robust at changing Rayleigh 
and/or position in the phase diagram (with more or less 
bubble nucleation in the system). This needs to be in- 
vestigated with new computational effort and will be re- 
ported in the future. 

Small-scales properties. Buoyant bubbles brings 
information from the physics of the bottom boundary 
layer in the bulk of the system. We then expect also 
in the bulk an increase of small scales fluctuations con- 
cerning both velocity and temperature flelds. In flg. (|3| 
we show the structure functions for vertical velocity and 
temperature: 



S^^FKt) = {[u,{^ + r) - ii,(x) • r]P),uik 

)bulk 



4p)(r) = ([r(x + r)-T(x)-fF) 



(10) 



where the average is restricted on all points x in the bulk 
of the cell and the increment r is always taken in hori- 
zontal directions. In the two panels of flg.(|3]) we show the 
results for both quantities for p = 2. For both flelds we 
have a viscous range very well resolved, where the struc- 
ture functions goes as oc r^. Therefore, the presence of 
large-bubbles do not destroy the differentiability at small 
scales, another signature that the numerical set-up is un- 
der control. Second, boiling system have an enhanced 
signal at small scales, meaning that energy dissipation is 
globally increased. Third, for the boiling case we start to 
see an inertial range with K41-like scaling oc r^/^. In the 
inset of both panels we measure the Flatness (or Kur- 



tosis) of each fleld i^n.,T(r) = sl^lj.{r) /[Sl^^J^j.{r)]^ at 
different scales, e.g. a way to quantify how much the 
pdf is close/different to a Gaussian. Intermittency, as 
measured by the deviation of the flatness from its Gaus- 
sian value, K = 3 becomes more and more important at 
decreasing the scale, in agreement with the general ob- 
servation that bubbles induce an increase of fluctuations 
in the system. For temperature (top panel flg. [3| the 
inertial range behavior is much more singular than the 
case for velocity, due to the enhancement of temperature 
jumps between inside and outside bubbles. At difference 
from the case for velocity, where large scale pdf is indis- 
tinguishable from a Gaussian {Ku^ 3 for r ~ L), here 
temperature is more sensitive to the presence of bubbles 
also at large scale. 

In conclusion, we have proposed and validated a novel 
LBM to attack multiphase flows with full consistent def- 
inition of heat exchange in the system (latent heat). We 
have applied this scheme to study convection under boil- 
ing condition and we have studied the effects of nucleat- 
ing large bubbles at the bottom boundary layer on both 
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FIG. 3: 2nd order Structure Function for velocity (bottom) 
and temperature (top) vs r (in lattice units) for boiling and 
non-boiling systems. Inset: values of Flatness (same sym- 
bols)) 



single point observable (temperature proflle and heat- 
flux) and two-point correlation functions (structure func- 
tions). The latter, allowed us to assess also the impor- 
tance of bubbles on small-scales velocity and temperature 
fluctuations, indicating an enhancement of the deviations 
from Gaussian statistics with respect to the non-boiling 
case. We acknowledge useful discussion with R. Verz- 
icco. The COST Action MP0806 is kindly acknowledged. 
Support within the DEISA Extreme Computing Initia- 
tive is acknowledged (FP6 project RI-031513 and FP7 
RI-222919). 
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